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Abstract 

A slender-jet model for the pinching of a liquid column is considered in the 
limit of vanishing viscosity. We find the model to develop a singularity in the 
gradients of the local radius and the velocity at a finite thread radius, so it 
does not describe breakup. However, the observed steepening of the profile 
corresponds to experiments and simulations with fluids at low viscosity. The 
singularity has similarity form, which we compute analytically. The result 
agrees well with numerical simulations of the model equations. 
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I. INTRODUCTION 



Considerable attention has been devoted recently to the breakup of an axisymmetric 
column of fluid The breakup is driven by surface tension forces, which reduce the 

surface area by contracting the fluid thread until its radius goes to zero at a point. As a 
typical example, we show a drop of water falling from a faucet in Fig.|^ 

Close to the point of breakup, the interface looks like a cone attached to a nearly flat 
interface. This is characteristic for low viscosity fluids, where viscosity is important only in a 
small spatial region around the point of breakup. Details of the initial conditions or external 
forces like gravity are believed to have little impact on the very localized behavior close to 
pinch-off. Indeed experiments with or without gravity and for a variety of nozzle diameters 
show very similar shapes 0]. This is because surface tension forces become very strong 
near pinching, and drive very small amounts of fluid. Thus the very rapid motion close to 
breakup is separated dynamically from the motion on the scale of the nozzle diameter both 
in space and time. A proper measure of length and time are the local scales 

4 = (^V)/7 , U = {i^'p')h^ (1) 

which only depend on the properties of the fluid. Here u is the kinematic viscosity, 7 the 
surface tension, and p is the density of the fluid. In the case that the minimum radius hmin 
of the fluid neck is much smaller than i^, a universal pinching solution has been observed 
IHJ^. As the radius of the fluid neck goes to zero, surface tension, viscous, and inertial forces 
are of the same order. 

For a low viscosity fluid like water, however, i,y is only 100 A, so this asymptotic solution is 
hardly of relevance experimentally. Accordingly, it would be extremely desirable to develop 
a similarity theory valid in the range hmin ^ ^u- Assuming that all flow features of a 
hypothetical Navier-Stokes solution are of the same order as hmm, it is tempting to look at 
solutions of the inviscid (Euler) equation for that regime. Unfortunately, the Euler equation 
is known to exhibit spurious blow-up of the local vorticity, even starting from smooth 
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initial data. This problem can be avoided by only considering the subclass of solutions 
which are irrotational, and are thus described by potential flow. Numerical simulations of 



inviscid, irrotational flow driven by surface tension were used by a number of authors PJTO 
to describe pinching. After the minimum radius has reached a value of about 1/20 of its 
initial value, all simulations show an overturning of the profile. This means that the neck 
radius h{z) is no longer a single- valued function of the position along the axis z. Recently ||5| 
it was shown that this overturning can be understood as the convergence onto a universal 
similarity solution of the inviscid, irrotational equations. Similar findings were reported 



by Chen and Steen |TT| for the related problem of a soap film which drives a flow in the 
surrounding air. Neglecting the inertia of the film, this corresponds to the motion of two 
fluids of equal density with surface tension between them. 

The fundamental question thus remains whether viscosity can really be neglected. A 
Navier-Stokes computation at low viscosities [1^ showed no signs of overturning, although 



the initial neck radius was A-lO'^i^,. Experimentally, there is also no sign of overturning 

even in water, an event that would show up in the side view as a flat wall perpendicular to 
the axis. Clearly, other and more detailed simulations of low-viscosity fluids are needed to 
confirm this, but there is reasonable concern that the assumption behind the inviscid-flow 
calculations are fundamentally flawed. 

In |jl3| we use a slender-jet model, originally developed in [0, to investigate the possible 



effect of a very small amount of viscosity on the pinching of a liquid thread. The model 
16|,[7[| can be shown to capture all the leading-order contributions to the jet dynamics below 



the viscous scale i^. In addition, it was demonstrated in that by including the full mean 
curvature as a higher-order contribution to the pressure, the model can reproduce experi- 
mental proflles even when hmm ^ As an example, in Fig.^a few proflles calculated from 
the slender jet model are superimposed on the experimental picture J^. Thus although the 
model is formally not valid in the presence of large gradients, it is interesting to understand 
the limit of small viscosity in this case since it turns out to be analytically tractable. 



In we show numerically that the slender-jet model for small viscosity develops very 



sharp gradients both in the local radius and the velocity field. This causes viscous effects to 
become important in a very small region where gradients are large even though the minimum 
radius hmin is still much larger than ii^. Clearly this is because the small-viscosity dynamics 
produces fiow features which are much smaller than hmin- 

In the present paper we show that the inviscid equations indeed have a singularity which 
leads to a blow-up of gradients in finite time. This has been done by explicitly constructing 
a similarity solution which solves the equations for u = 0. We extend the numerical code 



used previously in [|I^ to even smaller viscosities to show that in the limit v ^ oi small 
viscosities the slender-jet equations always select this inviscid singularity. The blow-up 
occurs while the minimum radius is still finite and an analytical description of the local 
profiles is given. In particular the inviscid equations are not able to describe breakup, 
although the full curvature term is kept. 

If one follows the above argument, inviscid scaling theories P, [I^ , p!5| may fail even for 
arbitrarily small i^. On the other hand, our analysis is based on a simplified one-dimensional 
description of a liquid thread, which represents a serious shortcoming compared to the 
three-dimensional axisymmetric representation of In particular, overturning cannot be 
described within the model. To settle the question whether inviscid scaling is consistent, a 
very careful comparison with both experiment [§,0 and full Navier-Stokes simulations 
is necessary. In particular more carefully resolved Navier-Stokes computations are highly 
desirable. 

In the next section we introduce the one-dimensional slender-jet model which forms the 
basis of our analytical description of the inviscid singularity. We then present conclusive 
numerical evidence for the existence of a inviscid singularity at a time tc- Derivatives of the 
surface profile and of the velocity diverge like a power law as function of tc — t. In the third 
section we present an analytical theory of the inviscid singularity. The resulting surface 
profiles agree well with numerical simulations. 
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II. MODEL AND SIMULATIONS 



The main assumption underlying the model of axisymmetric free-surface flow to be con- 
sidered here is that the fluid motion is directed mostly in the axial direction. This allows 



to set up an asymptotic expansion [16| in the thread radius, which at leading order gives 
equations for the radius h{z,t) of the thread and for the velocity v{z,t), which only depend 
on the axial coordinate z. In what follows we will deal with the model introduced in 

dth = —vhz — (2a) 
dtv = -vv,-p, + — , (2b) 

^ " /l(l + /l2)l/2 " (l + /i2)3/2' (2c) 

where the index refers to differentiation with respect to the variable. The fields h{z,t) and 
v{z,t) have been non-dimensionalized using some fixed length scale L of the problem (like 
the radius of a nozzle). The length L can be combined with surface tension 7 and density 
p to make up a time scale T = (pL^/'jY^'^ and a velocity scale U = L/T. Every quantity to 
follow will be non-dimensionalized using these units. 

Equation (^) expresses mass conservation for a radially uniform velocity field. Conser- 
vation of momentum (2b) not surprisingly has the form of Burgers' equation in the inviscid 



limit, driven by surface tension forces which are proportional to the mean curvature (cf. 
(|cD). Along the same lines of reasoning the equations (0) were guessed by Lee Jl^ for 



u = 0. By including the full mean curvature in (^ we have gone beyond the leading order 
asymptotics to exactly reproduce the static shape of a hanging drop suspended from an 
orifice. As an additional benefit, the most dangerous short-wavelength instabilities of the 
leading order model p = 1/h have been removed. In fact, the leading order model is elhptic 
for i> = ||TB|,[T^ and is thus ill-posed as an initial value problem. 

In [0,|T^ a finite difference scheme was developed capable of simulating (Q) at very low 
viscosities. To resolve the small-scale structures we are interested in, it is crucial to use an 
adaptive scheme, both in time and space. The minimum thread radius and the maximum 



gradient of h were taken as predictors where additional spatial resolution was necessary. 
Thus in a typical run the grid spacing at the position of the inviscid singularity was 6 
orders of magnitude smaller than at the boundary of the computational domain. Since 
the equations at low viscosity are very sensitive to noise, grids with smoothly varying grid 
spacings had to be used, where the spacing did not change by more than 1 % from one grid 
point to the next. With these precautions, no numerical damping or dissipation had to be 
used, except for the physical viscosity z/. 

Since we are interested in the limit of small viscosity, it would be tempting to put u = 
in (|2bD directly. However, we found that as soon as the motion is sufficiently nonlinear, 
our scheme developed instabihties on the scale of the grid, which caused the code to break 
down. Thus the inclusion of the full curvature term in is not enough to stabilize the 



numerical scheme. Similar short-wavelength instabilities have also been reported in ||20 
using a finite-element approach. On the other hand, exceedingly small amounts of viscosity 
are sufficient to stabilize the scheme, even though the viscous term is smaller than the others 
by several orders of magnitude throughout the domain. In the following, when speaking of a 
numerical solution of the inviscid equations, we will always refer to the limit of zero viscosity 
at constant time. 

We also experimented extensively with other means of regularizing the inviscid equations, 
for example by using numerical viscosities as in [0. In the upwind differencing scheme 
introduced in 0, the numerical viscosity is proportional to the grid spacing. The hope was 
to develop a scheme which automatically converges to the inviscid limit as one increases the 
resolution. Indeed, if the grid is coarse, a numerical viscosity was often sufficient to remove 
instabilities. But with improved resolution we always found the instabilities to return. Thus 
keeping a finite viscosity turns out to be the only reliable and at the same time the most 
physical way of dealing with the instabilities. These results indicate that the system (0) 
might be an ill-posed initial-value problem, in spite of the short-wavelength reguralization 
introduced by (P^. If on the other hand the problem is well-posed, and the instabilities are 



a problem of the numerical scheme, the above limit of small viscosity will yield a solution 
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which coincides with the one defined by the inviscid (Lee's) equations. 

It is clear that the slender-jet approximation represented by (0) assumes the existence of a 
single- valued function h{z) and is thus not able to describe overturning as observed for three- 
dimensional axisymmetric, irrotational flow P,p!0|,|5[] . Thus there is no direct correspondence 



to the pinching solution found in [Q. It could be that the inviscid singularity described in 
this paper is a consequence of the three-dimensional equations overturning. On the other 
hand, as explained in the Introduction, our working assumption is that the inviscid limit 
of the Navier-Stokes equation is singular, in which case the limit of vanishing viscosity has 
no connection with solutions of the Euler equation. Thus the fact that full the potential 
flow equations show overturning has no obvious implications for the limit of of vanishing 
viscosity studied here. The resolution of these issues must clearly await full Navier-Stokes 
computations at small viscosity. 

Figure ^ shows a simulation of (0) at a very large Reynolds number in a liquid bridge 
geometry with initial radius r = 0.05 and a small sinusoidal perturbation of wavelength 
A = Atttq superimposed on it. Owing to the Rayleigh instability, the bridge starts to pinch. 
Shown are three proflles close to the inviscid singularity, where the minimum radius has 
already decreased by a factor of 20. To the right of the minimum, an almost conical neck 
region is seen, on the other side a round drop has formed. Because of its small radius, the 
pressure in the flat region is high, pushing fluid over to the right. This causes the interface 
to perform a sliding motion, which lets the interface become steeper and steeper, since the 
drop cannot move appreciably owing to its large inertia. Note that the pressure goes to a 
value close to zero at the eventual place of the inviscid singularity, marked by an arrow. 

Next we zoom in on the point around which hz goes to inflnity, marked by an arrow in 
Fig.0. To sufficiently resolve the singularity without letting viscosity become important, the 
Reynolds number is chosen to be Re = 4.7 • 10^, much larger than the Reynolds numbers 



treated in [|T^. With such a small Reynolds number we are no longer able to resolve the 
huge range of length scales between the outer and the viscous scale, because L/i^ is now 
2 • 10^^. However, the early stages of the evolution of the liquid bridge, where viscous effects 



are still small can safely be resolved. In |T3| it was demonstrated that the slope goes to 



infinity near the inviscid singularity and thus the singularity time tc can be estimated from 
the blow-up of hz. It follows from our scaling theory, to be presented in the next section, 
that {hz)max ~ {ic — t)"^ ■ Thus tc can be computed very accurately by plotting {{hz)max)~^ 
versus time and fitting with a linear law. In Fig.^ we plot the maximum pressure gradient 
Pz) which drives the fiuid motion, and the maximum velocity gradient as a function of tc — t. 
It is seen that both pz and Vz settle on a power law 

{Pz)max ~ {ic — {Vz)max ~ {tc — t)~^ . (3) 

The pressure gradient contains the highest (third) derivative in the problem. The fiuctu- 
ations seen in the curve thus give an estimate of the amount of noise introduced by the 
regridding procedure. No noise is seen in the first velocity derivative, which clearly confirms 
the scaling given in (0). We thus see that the inviscid singularity is governed by power law 
scaling, which will be investigated in more detail in the next section. We made sure that for 
all of the evolution shown the viscous term remained small compared to other terms. Thus 
the observed scaling represents the inviscid limit of the equations we want to investigate. 

We also considered a great variety of other initial perturbations, which changes the initial 
evolution of the bridge. In addition, we considered the bifurcation of a hanging drop in the 
presence of gravity as shown in Fig. ^ In all these different cases, we always found the same 
singularity with scaling (j^) to be selected in the inviscid limit. Thus there is good evidence 
that the inviscid singularity to be described below is a universal feature of the equations, 
independent of initial conditions. 

The equations (0) with {u = 0) or systems very similar in structure have been used by a 
number of researchers [ p^ , p!^ ,pl|- p3| , |20| , p!8| , |2^ to describe pinching. However in Lee's original 



paper and most of the later work no attempt is made to resolve any detailed structure in 
the pinch region. For example, in units of the length of the computational domain, the grid 
spacing is (ix = 1/20 in |]T^ and 1/50 in For comparison, the minimum grid spacing 



used in the present paper is dxmin = 10 In most papers [Il7| , ^ , p2| , p0| the computation 
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is stopped at hmin between 1/10 and 1/20 of the initial radius, when the position of the 



minimum and the drop size can be faithfully estimated. In ||2^ , which uses a finite element 
code, the computations were stopped when numerical instabilities on the scale of the grid 



were observed. A spectral method is used in [Q, but the computation was stopped at 
hmin = 0.3 due to numerical problems. Thus our finding of a short-wavelength instability of 
the inviscid equations is consistent with other numerical work, using a variety of different 



methods. It remains to be seen whether the inviscid versions of the Cosserat model [21 



or 



other asymptotic models 24] are more well behaved. 



III. INVISCID SIMILARITY SOLUTION 

We have seen in the previous section that derivatives grow sharply near the inviscid 
singularity, while on the other hand the height and the velocity remains finite. This means 
that the similarity ansatz of [jT| has to be generalized to include a "background" height and 
velocity profile, which is slowly varying on the scale of the singular part. At the same time, 
the singularity may be moving with some speed Vg, which is not necessarily the speed V 
with which it is convected. Thus one ends up with the similarity form 
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h{z',t') = H + rf('''^^'^'^ 



v{z',t') = V + t"^g 



(4) 

'z' + Vst'\ 



where 

z = z — Zc and t' = tc — t (5) 

measures the spatial and temporal distance from the singularity, respectively. On the spatial 
scale on which the inviscid singularity develops, H, V, and Vg are approximately constant. 
Note that @ has self-similar form, which is superimposed on a traveling wave solution. 
Also, we assumed h{z',t') and v{z',t') to have the same scaling exponents, because this 
automatically balances the terms hzV and Vzh in (pa]) . 

For the ansatz to be consistent, one needs a > 0, so in the singular limit t' — > one 
is left with the finite height H and velocity V. For the derivatives to blow up, (3 > a. To 
see whether (|D solves the model equations (|^), we balance the most singular terms in t', 
deriving equations in the similarity variable 

z' + Vst' 



One thus finds from 



iPf'v-af)t 



/a-l 



(V - Vs)f + \Hg' 



t""-". (7) 



We expect the right hand side of to make the dominant contribution, which will be 
the case if (3 > 1. This is because then the function / drops out of the equation, and 
the similarity equations only depend on the derivatives /' and g'. Thus both / and g are 
determined only up to constants, which is needed for consistency because our ansatz (Q) has 
a free constant built in. 

Consequently, the angular bracket must vanish, giving 

9 = jj / , (8) 
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which means that up to constants and a difference in amphtude 2{Vs ~V)/H the profiles of 
the height and of the velocity are the same. 

Turning to (pE|) with z/ = 0, one finds to leading order 



- (K - V)g't'^-^ = 4- t'-'". (9) 



' ~[fV 

Note that the term on the left corresponds to the highest derivatives in as given by 
Balancing the left and the right hand side, one finds the scaling law 

p = 3a. (10) 

Combining (H) and the similarity equation reads 

-a/'=(f)', (11) 
where a = 2{Vs — VY/H. Evidently, the constant a can be eliminated by the transformation 

0(r^) = a-VY(r/). (12) 
The most general solution of the equation for 0, 

has the form 

4>{r,) = ct>,F[ct>T{r]-Vo)i (14) 



where F[^) is a particular solution of 

Equation (0) can easily be solved using standard tricks. In view of the freedom 



implied by (|T4D, we choose -F(0 to have its maximum at ^ = 0, and to fall off to 1/2 at 
^ = ±1/2. Then F is given implicitly by 

e = ^(l + 2F)(l-F)V2. (15) 

This function is represented in Fig. ^ as the solid line. It decays to zero like ^^^''^ at infinity. 
In view of the similarity form (^ this ensures that the leading dependence on t' drops out 
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far away from the singular point. This is necessary for the solution to match onto the slowly 
varying background field. 

We note that the singularity described above is not just kinematic in nature, since a 
contribution from the capillary forcing enters the dominant balance in (|^). In addition, the 
form of the singularity with K 7^ ^ implies that it is also a traveling wave. Surprisingly, 
the local shape (|l^) of the first derivative of the local radius is identical to that of a scaling 
solution h = t'G{z' /t''^/'^) of the simple kinematic wave dth + hdzh = This can be 
checked directly from the similarity equation for G{^). At present we do not know if this is 
a coincidence or is the result of a deeper analogy. 

To test the prediction of the theory, we took the same simulation as in Fig. ^ at a 
time where the slope was = 10^. We included the profile of in Fig. ^, shifting the 
maximum to the origin and normalizing its width. It is evident that our analytical theory 
is fully confirmed by the comparison. The free constants (po, rjQ in (|1^ and the parameter a 
in ( p!TD are not determined by the similarity theory. Indeed, we confirmed that they depend 
on initial conditions and therefore cannot come out of a local theory. 

Our next concern is to find the scaling exponents a and /3, which are not determined 
from dimensional reasoning as in [0], but which are constraint by the scaling relation (|rU|). 
In addition, ([Tl|) does not depend on the values of the exponents, so a cannot be selected 
by properties of the similarity equation, as was the case in [^. To investigate this prob- 
lem, we must look at next to leading order terms such as the ones contained in Eq. (^. 
Correspondingly, there are sub- leading terms in h and f , which have the form 

h{z\ t')=H + n + A + . . . , (16) 

and correspondingly for v{z',t'). Then (pa]) becomes to next to leading order: 



[3af'v - af] t"^-' = -gf - \fg' + (K - V)f[ - 



±l—a 



(17) 



Since the terms must balance, we get 
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which are the desired exponents. Note that this conforms with the scahng of both the 
pressure and the velocity gradient from since 

P. ~ t'-'" (^1^) and Vz^t"'-^g'. (19) 

Thus both the value of the exponents and the shape of the profiles is in excellent agreement 
with theory. 



IV. DISCUSSION 

In we have studied the steepening of the height profile for small but finite viscosities. 
The slopes saturate at a large value, the maximum slope roughly following a scaling law 
{hz)max ~ Re^'^^. As long as the viscous term is much smaller than the pressure gradient, 
one can use the inviscid similarity solution. The naive expectation is that the slope saturates 
as soon as the viscous term is of the same order as the pressure gradient. Using (^) and 
(p!8|), the temporal scaling of the pressure gradient is Pz ~ t''^ and that of the viscous 
term Re~^ {v zh) z / h"^ ~ i?e~^t'~^/^, where Re~^ = u/{UL) is constant. Equating the two we 
find {hz)max ~ -Re^/^, which is far too small an exponent. A possible explanation is that the 
Reynolds numbers for which the exponent 1.25 was found are still too small. But more likely 
there is an intricate interplay between the inviscid singularity and viscosity, leading to a more 
complicated intermediate scaling range. Indeed, for the Reynolds numbers considered, the 
slope continues to grow far beyond the time where the pressure gradient first balances the 
viscous term at a point. It thus remains a challenge to find the mechanism which makes the 
slope saturate. 

The pinching of vortex sheets and jets within the framework of inviscid hydrodynamics 
has recently been the object of much scrutiny p8|-p0|]. The motion is either driven by 
surface tension, or comes from the flow field generated by the regions of high vorticity. Short 
wavelength instabilities arise just as in our problem, which are usually removed by filtering 
out high-wavenumber components of the Fourier spectrum It would be interesting to 



13 



see whether this is equivalent to adding a small amount of viscosity, which is the physical 
regularization used in the present paper. It is not clear that Krasny filtering actually selects 
the physical solution. So far, analytical solutions have not been found for the singular motion 
of vortex sheets or lines. It remains to be seen whether the singularity structure found in 
the present problem carries over. 

We have seen that the slender-jet model (0) in the limit of small viscosity is characterized 
by by more than just the scale of the minimum radius. Instead, a shock- type singularity de- 
velops whose width represents another, much smaller scale. The crucial question is of course 
whether a similar mechanism is at work in the small-viscosity limit of the Navier-Stokes 
equation, which could make this limit singular. It is unlikely that the three-dimensional 
equation has precisely the same spatial singularity structure than the one found in the 
model equations, which constrain axial velocity gradients to a far greater extend. Instead it 
is probable that high-pressure fluid in the neck is injected into the drop, a situation which is 
only poorly captured by the slice average of the model. The Navier-Stokes equation might 
thus form a very thin boundary layer around the jet, in which viscosity remains important 
even for arbitrarily small v. This is a possible scenario which would invalidate the purely 
inviscid calculation of 0, at least form a physical point of view. 

Even for the model equations, many unanswered questions remain. It is fascinating 
that even a simple one-dimensional model is capable of such complexity, the key to its 
understanding lying in the analysis of the singularities. 
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FIGURES 



FIG. 1. A drop of water falling from a faucet 0.51cm in diameter. The lines represent a 
computation using one-dimensional model equations. This figure is taken from Q. 

FIG. 2. A closeup view of simultaneous radius, velocity, and pressure profiles close to the 
inviscid singularity. The Reynolds number is Re = 5 • 10^, and the time between successive profiles 
IS 5 • 10"^ The arrow marks the asymptotic location of the shock. 

FIG. 3. The maximum of the gradient of the pressure and of the velocity as function of 
the time distance from the inviscid singularity. The Reynolds number is 4.7 • 10^. Both curves 
asymptote to a slope very close to —1. 

FIG. 4. The normalized similarity function F{^), cf. (|l^), ([l^). 
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Figure 1 
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Figure 3 
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